#ifndef ArrayBasedAxisAlignedDistribution_hpp
#define ArrayBasedAxisAlignedDistribution_hpp

#include "functors/IFunctor.hpp"
#include "domains/IAxisAlignedDomain.hpp"
#include "distributions/IAxisAlignedDistribution.hpp"
#include "linalg/Vector.hpp"
#include "other/utils.hpp"		// for binary interval search

template <class UFUNC, class VFUNC>	// the template is forced here since we can't keep as members and subsequently call operator() on the IFunctors needed to specify the 2D distribution function
class ArrayBasedAxisAlignedDistribution : public IAxisAlignedDistribution {	// UFUNC and VFUNC must derive from IIntegrableFunctor
	public:
		ArrayBasedAxisAlignedDistribution();
		~ArrayBasedAxisAlignedDistribution();
		virtual bool sample(double &_u, double &_v, double cdfvalu, double cdfvalv);
		virtual void setBounds(double umin, double umax, double vmin, double vmax);
		void setResolution(unsigned _nptsu, unsigned nptsv);
		void setFunctors(UFUNC _ufunc, VFUNC _vfunc);
		void build();
		//
		double get_axis1_min() { return umin; }
		double get_axis1_max() { return umax; }
		double get_axis2_min() { return vmin; }
		double get_axis2_max() { return vmax; }
	protected:
		IAxisAlignedDomain *dom;
		UFUNC ufunc;
		VFUNC vfunc;
		double umin,umax,vmin,vmax;
		unsigned nptsu,nptsv;
		double *u,*v,*cdfu,*cdfv;
};

template <class UFUNC, class VFUNC>
ArrayBasedAxisAlignedDistribution<UFUNC,VFUNC>::ArrayBasedAxisAlignedDistribution() {
	u=NULL;
	v=NULL;
	cdfu=NULL;
	cdfv=NULL;
	nptsu=0;
	nptsv=0;
	dom=NULL;
}

template <class UFUNC, class VFUNC>
ArrayBasedAxisAlignedDistribution<UFUNC,VFUNC>::~ArrayBasedAxisAlignedDistribution() {
	if (nptsu>0) {
		delete [] u;
		delete [] cdfu;
	}
	if (nptsv>0) {
		delete [] v;
		delete [] cdfv;
	}
}

template <class UFUNC, class VFUNC>
void ArrayBasedAxisAlignedDistribution<UFUNC,VFUNC>::setBounds(double _umin, double _umax, double _vmin, double _vmax) {
	umin=_umin;
	umax=_umax;
	vmin=_vmin;
	vmax=_vmax;
}

template <class UFUNC, class VFUNC>
void ArrayBasedAxisAlignedDistribution<UFUNC,VFUNC>::setResolution(unsigned _nptsu, unsigned _nptsv) {
	nptsu=_nptsu;
	nptsv=_nptsv;
}

template <class UFUNC, class VFUNC>
void ArrayBasedAxisAlignedDistribution<UFUNC,VFUNC>::setFunctors(UFUNC _ufunc, VFUNC _vfunc) {
	ufunc=_ufunc;
	vfunc=_vfunc;
}

template <class UFUNC, class VFUNC>
void ArrayBasedAxisAlignedDistribution<UFUNC,VFUNC>::build() {	// TODO: somewhere between build() and sample() we lose resolution in u- and v-space using MaxwellianFunctors with vth<<vd (probably not a big issue for most plasma sims)
	u=linspace(umin,umax,nptsu);
	v=linspace(vmin,vmax,nptsv);
	double *fu=new double[nptsu];
	double *fv=new double[nptsv];
	cdfu=new double[nptsu];
	cdfv=new double[nptsv];
	//
	for (unsigned i=0; i<nptsu; i++) {
		fu[i]=ufunc(u[i]);					// assume here that f[i]>0 always, and ignore normalization
		if (i==0)
			cdfu[i]=0;
		else
			cdfu[i]=cdfu[i-1]+(fu[i]+fu[i-1])/2.;	// using trapezoidal rule
//		cdfu[i]=ufunc.integrate(u[0],u[i]);
	}
	double dcdfu=cdfu[nptsu-1]-cdfu[0];
	for (unsigned i=0; i<nptsv; i++)
		cdfu[i]=(cdfu[i]-cdfu[0])/dcdfu;		// map cdf to interval [0,1]
	//
	for (unsigned i=0; i<nptsv; i++) {
		fv[i]=vfunc(v[i]);					// assume here that f[i]>0 always, and ignore normalization
		if (i==0)
			cdfv[i]=0;
		else
			cdfv[i]=cdfv[i-1]+(fv[i]+fv[i-1])/2.;	// using trapezoidal rule
//		cdfv[i]=vfunc.integrate(v[0],v[i]);
	}
	double dcdfv=cdfv[nptsv-1]-cdfv[0];
	for (unsigned i=0; i<nptsv; i++)
		cdfv[i]=(cdfv[i]-cdfv[0])/dcdfv;		// map cdf to interval [0,1]
	//
	delete [] fu;
	delete [] fv;
}

template <class UFUNC, class VFUNC>
bool ArrayBasedAxisAlignedDistribution<UFUNC,VFUNC>::sample(double &_u, double &_v, double cdfvalu, double cdfvalv) {
	unsigned idxu,idxv;
	bool foundu=binaryIntervalSearch(cdfu,nptsu,cdfvalu,idxu);
	bool foundv=binaryIntervalSearch(cdfv,nptsv,cdfvalv,idxv);
	if (foundu && foundv) {
		if (idxu<nptsu-1)
			_u=u[idxu]+(cdfvalu-cdfu[idxu])/(cdfu[idxu+1]-cdfu[idxu])*(u[idxu+1]-u[idxu]);
		else
			_u=u[nptsu-1];
		//
		if (idxv<nptsv-1)
			_v=v[idxv]+(cdfvalv-cdfv[idxv])/(cdfv[idxv+1]-cdfv[idxv])*(v[idxv+1]-v[idxv]);
		else
			_v=v[nptsv-1];
		//
	}
	return false;
}

#endif
